Sampling-Based Resolution-Complete Algorithms 



In this paper, we describe a novel approach for checking safety speci- 
fications of a dynamical system with exogenous inputs over infinite time 
horizon that is guaranteed to terminate in finite time with a conclusive 
answer. We introduce the notion of resolution completeness for analysis 
of safety falsification algorithms and propose sampling-based resolution- 
complete algorithms for safety falsification of linear time-invariant discrete 
time systems over infinite time horizon. The algorithms are based on de- 
terministic incremental search procedures, exploring the reachable set for 
feasible counter examples to safety at increasing resolution levels of the 
input. Given a target resolution of inputs, the algorithms are guaranteed 
to terminate either with a reachable state that violates the safety specifi- 
cation, or prove that no input exists at the given resolution that violates 
the specification. 

1 Introduction 

1.1 Background 

Simulation-based techniques for formally verifying properties (called specifica- 
tions) of discrete, continuous and hybrid systems, have come under great deal 
of attention recently. The motivation to use such techniques arises from the 
fact that most of the real-world systems are quite complex and operate in the 
presence of unknown external disturbances. As a result, verifying that each 
and every state of a system satisfies a given specification may be impractical, 
or in general even impossible. The problem of finding the set of all states the 
system can reach (called as the reachable set), based on its dynamics and initial 
conditions, is known as the reachability problem in literature. For continuous 
and hybrid systems, this problem is in general known to be undecidable [1, 2]. 
An important class of specifications are safety specifications that describe the 
properties that the state of a system should satisfy, to be considered safe. For 
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analyzing safety specifications of a system, a wide variety of methods have been 
proposed in [3, 4, 5, 6, 7, 8, 9, 10, 11]. 

Most of these methods attempt to verify safety of a given system by over 
approximating the actual reachable set. Hence, they are liable to generate a 
spurious counter example which violates a specification, but is not a feasible 
trajectory [12]. Even though refining the abstraction is usually possible, there 
is in general no guarantee that the process of successive refinements would stop 
in finite time [12]. As a result, such methods can only verify safety of a system 
but will be inconclusive with regard to disproving it. Safety of a given system 
can be disproved only by working with either the actual reachable set, or, by 
giving a feasible counter example (constructed, e.g., using a simulation-based 
method) . 

1.2 Sampling-based algorithms for safety falsification 

To answer the complementary question of safety falsification, sampling-based in- 
cremental search algorithms have been proposed by us, and others, in [13, 14, 15, 
16], that are based on similar algorithms used in robotics [17]. These algorithms 
try to falsify safety of the system quickly, but they are only probabilistically- 
complete (see [17]). This means that, the algorithms will find a counter example 
(if one exists) with probability 1, if they are allowed to run forever. However, 
if they are terminated before a counter example is found, then they become 
inconclusive. One such instance is shown in Fig. 1, at the left. The trajectories 
are generated using Rapidly Exploring Random Tree [17], a probabilistically- 
complete algorithm, for a point mass moving with bounded velocity in two 
dimensions starting from a set containing origin. Terminating the search pro- 
cedure before the unsafe set is reached leads to the wrong conclusion that the 
system is safe. Note that the samples are points (zero volume sets) and the 
sampling procedure is randomized. 

To analyze the completeness properties of search based motion planning algo- 
rithms used in robotics, the notion of resolution completeness has been proposed 
in [18, 19]. A resolution-complete algorithm, is guaranteed to find a solution 
(if one exists), in finite time, provided that, the resolution of discretization in 
input space, and state space, is high enough. In [20, 21], we introduced a similar 
notion for disproving safety of continuous and hybrid systems over infinite time 
horizon, and, proposed resolution-complete deterministic algorithms applicable 



to linear time invariant systems (abbreviated as LTI systems). The algorithms 
work by incrementally building trajectories in state space at increasing levels 
of resolution, such that, cither they fetch a legitimate counter example, or a 
guarantee that no such counter example exists at given level of resolution (and 
hence a conclusive answer to the unsafety problem), in finite time. In Fig. 1, at 
the right, we show an example where we end up with a non-zero volume under- 
approximation to reachable set (when no counter example was found), when 
using a sampling-based resolution-complete algorithm for safety falsification of 
a second order system with exogenous inputs [20] .Very recently, an alternate 
notion of resolution completeness has also been proposed by Cheng and Kumar 
for continuous-time systems for the case of finite horizon in [22] . 

1.3 Contributions of the paper and relation to other ap- 
proaches 

In this paper, we propose two new resolution-complete algorithms that use in- 
cremental grid-based sampling methods (similar to those proposed in [23]) with 
good coverage properties for exploring the state space. The first algorithm uses 
breadth-first-search based scheme with branch and bound strategy to explore 
the state space for counter examples to given safety specification, at increasing 
levels of resolution of the input. This algorithm can be applied to discrete- 
time LTI hybrid systems. Simulation results indicate that this algorithm, is an 
improvement over the one proposed by us in [21] which is based on depth-first- 
search based scheme with branch and bound strategy. The second algorithm 
proposed in this paper can be used to falsify safety of discrete-time LTI con- 
tinuous systems more efficiently, when the initial set is the equilibrium point. 
The reachable set for such initial conditions is a convex set. The proposed algo- 
rithm uses this fact to explore the state space more efficiently. Since both the 
algorithms are resolution-complete, they consider the safety falsification prob- 
lem over infinite time horizon, with guarantees of finite-time termination of the 
search procedure and a conclusive answer at termination. 

An important feature that distinguishes our approach from alternate ap- 
proaches recently proposed by others (e.g. [24]) for addressing safety over infi- 
nite time horizon, is that the requirements on discretization of state space in our 
case do not depend on time length of trajectories. This is an important advan- 
tage by itself, and also helps the algorithms in avoiding the so called wrapping 
effect [25]. This means that, in case no counter example is found, then the 
quality of the approximation constructed as a proof for safety of the system is 
not affected by time horizon. We also do not discretize the space of inputs to 
obtain completeness guarantees (unlike the approaches presented in [22, 24]). 

The paper is organized as follows. We introduce the framework for describing 
hybrid systems and reachable sets in Section 2 together with a formal definition 
of the notion of resolution completeness for safety falsification. In Section 3, we 
explain the main idea used in the proposed algorithms for resolution complete- 
ness. Conditions for resolution completeness of the proposed algorithms are 
discussed in Section 4. In Section 5, we explain the sampling method used by 



us in the algorithms and in Section 6, we present the algorithms, together with 
the proof of their completeness. Simulation results are discussed in Section 7, 
and the paper is concluded in Section 8. 

2 Preliminaries 

In this section, we introduce notation for describing hybrid systems and reach- 
able sets and define the notion of resolution completeness. 

Definition 1 (Hybrid System). We define a discrete-time LTI hybrid system 
H as a tuple, H = (Q, X,U, U, <&, A,X, S,T), where: 

• Q is the discrete state space. 

• X C MP is the continuous state space. 

• U is the family of allowed control functions equipped with a well defined 
metric. Each control function u € U is a function u : [0,tf] — > U, where 
tf 6 N is the terminal time of the trajectory. The convex set U C K m is 
the input space. For simplicity, we assume U to be a unit hypercube. 

• &:QxXxU—>Xisa function describing the evolution of the system on 
continuous space, governed by a difference equation of the form x(i + 1) = 
3>(x, q,u) — A q x{i) + Bqu(i),i € N, and A qi B q are real matrices of size 
n x n, n x m respectively. 

• Ac (Q x X) x (Q x X), a relation describing discrete transitions in 
the hybrid states. Discrete transitions can occur on location- specific sub- 
sets Q(q,q') C X, called guards, and result in jump relations of the form 
(q,x) h-> (q',x). 

• I,5,T C Q x X are, respectively, the invariant set, the initial set, and 
the unsafe set. 

The semantics of our model are defined as follows. When the discrete state 
is in location q, the continuous state evolves according to the difference equation 
x(i + 1) = A q x(i) + B q u(i),i € N, for some value of the input u(i) € U, with 
(q(0),x(0)) € S. In addition, whenever x(i) G Q(q(i),q'), for some q' , the 
system has the option to perform one of the discrete transitions modeled by 
the relation A, and be instantaneously reset to the new discrete state q' , while 
the continuous state remains the same as before the discrete transition. The 
system is required to respect the invariants by staying within X at all times. 
For the cases when the discrete state space is just a single location, we drop 
the discrete state q from the notation, ft — {q, q : [Q,tf] — > Q} denotes the set 
of trajectories on the discrete space. Trajectories of the system starting from 
z € Q x X and using u € U, under the discrete evolution q, are denoted by 
^(z, u, q) C Q x X. A point on the trajectory ^(z, u, q), reached at time i < tf, 
is denoted by ip(z,u,q,i) € Q x X. We will denote by A° the interior of a set 
A. 



Definition 2 (Reachable Set). The reachable set 1Z(U) for a system H denotes 
the set of states that can be reached in the future. It is defined as 

1Z{U) = |J *{z,u,q). 
zes, u&j, gen 

We now formalize the notion of resolution completeness of an algorithm for 
safety falsification of a discrete-time LTI system with exogenous inputs. 1 

Definition 3 (Resolution completeness). A given algorithm is resolution-complete 
for safety falsification of a system H , if there exists a sequence of family of con- 
trol functions, {UjYjLx, satisfying Uj cUj+i,Vj, and liuij^ooUj = U (in the 
sense of a given metric), such that, for any given j > 1, the algorithm termi- 
nates in finite time, producing, either a counter example ip(za,u,q,t), using a 
control function u G 14, zq G S, q G fi, t > 0, satisfying, ip(zo,u,q,t) 6 T, or a 
guarantee that, TZ°(Uj) n T = 0. 

3 Basic Idea 

One way to achieve completeness is to construct an approximation IZj (while 
searching for a counter example) that satisfies the set inclusion lZ°(l4j) C IZj C 
1Z(U) (shown in Fig. 3(a)), thus guaranteeing feasibility of counter examples and 
safety with respect to control functions belonging to 14 j . The algorithms that we 
propose in this paper use this idea. To construct IZj, the algorithms discretize 
the state-space using multi-resolution grids. For a discrete location q G Q, 
G(q) denotes the multi-resolution grid for location q that over-approximates 
1Z n T(q, •). The algorithms keep a record of the portion of G(q) that is found 
to be reachable (denoted as Gf(q)), either a priori or during an execution of 
the algorithm. G u (q) denotes the rest of G(q), i.e. G u (q) = G(q) \ Gf(q). The 
algorithms progressively sample regions of G u (q) (called cells) at increasing 
levels of resolution and try to construct trajectories that end in the sampled cell 
starting from somewhere in Gf(q). The conditions for state-space discretization 
derived in Section 4 guarantee that finding one feasible point ip(zQ,u,q,t) in a 
cell €(ej(q)) of size Sj(q), with z G Gf(q),u G Uj,q G ft, is enough to claim 
that £ C 1Z(U). In Fig. 2, we show an execution of such a resolution-complete 
safety falsification algorithm (starting at resolution j and stopping at j + 1) for 
the case when card Q = 1 . 

We would like to remark here that, if we can find conditions that ensure 
the set inclusion, TZ°(Uj) C IZj C 1Z(U) for nonlinear systems, then similar 
algorithms can be used for resolution complete safety falsification of non linear 
systems as well. 



1 For a similar notion applicable to more general class of systems, we refer the reader to [21]. 
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Figure 2: Execution of a resolution-complete algorithm for card (Q) = 1 

4 Conditions for Resolution-Complete Safety Fal- 
sification 

In this section, we derive necessary conditions for resolution-complete safety 
falsification of discrete-time LTI hybrid systems. For all the discussion that 
follows, ||-|| denotes the infinity norm. The set of control functions is, U — 
{u{-) : \\u{i)\\ < 1, Vi G N, i < t f }. C q = [B q A q B q . . . A 1 q l ~ 1 B q ] denotes the 
controllability matrix of the system in location q G Q. 

4.1 Control functions for resolution completeness 

For resolution completeness 2 , we need a sequence of control functions 
such that Uj C Uj+i,Vj G N and lim^oo Uj = U in some metric defined on 
the space of control functions. We consider sequence of families of piece- wise 
constant control functions for our algorithm. 

Proposition 1. For given family of control functions U , the sequence of family 
of control functions {Uj}^, Uj = {«(•) : \\u(i)\\ < lj, Vi G N, i < t f ,t f G N} ; 
with {lj} a strictly non decreasing sequence of real numbers and lim^oo lj = 1 
satisfies Uj C Wj+i, Vj and lim^oo — U in norm. 

Proof. The proposition is proved by stated requirements on {lj}. M 

4.2 Assumptions 

Assumption 1. For each discrete location q G Q, the system H is stable at the 
origin, rank (B q ) = m, and, rank (C q ) = n. 

Stability (along with Assumption 2) guarantees that G(q) has a finite volume, 
rank (B q ) = m and rank (C q ) = n is needed to be able to use Propositions 2, 3. 

Assumption 2. For each discrete location q G Q, S(q, •), G{q, •) are specified 
as convex polytopes andX(q, -),T(q, •) are specified as a convex polyhedra. 

2 Resolution is defined on the space of control functions, Li. 
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Figure 3: Set inclusions for resolution completeness 

At each step, the algorithms incrementally build the trajectories and check 
for unsafety and discrete mode switches by solving linear programs. Hence we 
need the sets S{q-),l{q, -),G{q, -),T{q, •) to be convex polyhedra. Boundedness 
of S{q, ■), G{q, ■) is used to prove finite time termination of the algorithms. 

We now consider continuous dynamics in a given discrete location 5 £ Q, 
and derive sufficient conditions for discretization of cells in G{q). 

4.3 Sufficient conditions for state space discretization 

We first state an important proposition that guarantees that the set of points 
ip{z, u, q, k), that can be searched by the algorithms for feasibility, starting from 
z, with u £ Uj (as in Proposition 1) in k £ N steps has a non zero volume. This 
follows from the assumption that rank {B q ) — m, and, rank (C q ) — n in each 
discrete location q £ Q. 

Proposition 2. For a given system H satisfying Assumption 1, and for a given 
location q £ Q, the set of points reachable by using a control function u £ Uj 
(as in Proposition 1) over k steps has a non empty interior if \njm\ < k < n 
for any j > 1, where n is the dimension of the continuous state space and m is 
the dimension of input space. 

Proof. Let k be chosen such that \n/m~\ < k < n. The continuous dynamics 
of the system are x{i + 1) = A q x{i) + B q u, where u £ Uj. Applying this over 
k steps, we get x(k + i) = A q x(i) + [B q A q B q . . . A^B^u, where u £ R km is 
the augmented input over k steps. Let v = [B q A q B q . . . A^ 1 B q ]u. Then the 
dynamics can be written as x(k + i) = A q x(i) + v,v£ R n . Since rank (C q ) = n, 
rank [B q A q B q . . . Az~~ 1 B q ] — n. Hence the set of points reachable in k steps is 
guaranteed to have a non empty interior. ■ 

We now derive a conservative upper bound on the discretization of G(q), so 
that for a cell £, if £ n 11° {Uj) nl(q, 0^0 then £ C 1Z{U) nl{q, ■). To do so, we 
first prove the result for a simpler case in Lemma 1, and then prove the main 
result in Propositon 3. 

Lemma 1. Consider a continuous system H , with x{i + 1) = Ax{i) + v with, 
x, v £ W l . For a given X\, and an ot\ > 0, with x\ — Axq + vi,Xq,Vi £ 
K n , and, \ \ vi\\ < ol\, the following holds true: For any X2 £ R n , and a given 



a 2 > &i, if \\xi — X2W < s, and e < a 2 — ot\, then 3t> 2 G R", such that, 
x 2 = Axq + v 2 , with, \\v 2 \\ < a 2 - 

Proof. Please refer to Fig. 3(b). V2 = v\ + X2 — x\ proves the result. ■ 

Proposition 3. Consider a system H satisfying Assumption 1, with sequence 
of families of control functions {Uj}JL 1 , as in Proposition 2. Let j e N and 
q G Q be fixed. Then, if the cell size £j(q) for a cell £ £ G(q) satisfies the bound, 
Ej(q) < (l-lj)/\\T+\\, whereT q = [B q A q B q . . . A q ~ 1 B q ] andT q + is the pseudo 
inverse, then £ nlZ°(Uj) n% ■) ^ =*> £ C 1l(U) n% •)■ 

Proof. The dynamics of the system over k steps can be written as x{k + i) = 
A k q x{i) + v, with v = T q u. T q = [B q A q B q . . . A q ~ 1 B q ] and u € R km , \\u\\ < I 
is the augmented input over k steps. This implies that u — T^v, where is 
the pseudo inverse of T q . Finding the tightest bounds on v is hard. However 
a conservative bound on v is ||u|| < l/T + . Now, let vi = T q iii and w 2 = r 9 7i 2 
with ui,u 2 € M. km and ||ui|| < lj and ||u 2 || < 1- Let lZ(z, k) denote set of points 
reachable by the system by using u £U with tj = k, and T^°{z, k) the interior 
of the set of points reachable by the system by using v! € Uj with tf = k, under 
continuous evolution, starting from z. This is shown in Fig. 3(c). The result of 
Lemma 1 implies that for e q < (1 - l 3 )/T+, £ n U°(z,k) ^ => £ C K(z,k). 
This proves that £ n TZ^Uj) n l(q, ■) ^ ^ £ C ^(W) n T(g, •)■ ■ 

5 Incremental Grid Sampling Methods 

As we discussed in Section 3, for each location q G Q, we use a multi-resolution 
grid. Assume for this section that, card (Q) = 1, and that, G is a rt-dimensional 
unit cube, whose origin is the origin of coordinate axis. Using an iterative ly 
refined multi- resolution grid ensures that for a given j, G u with resolution Sj 
does not have to be built from scratch. For our work, we will use multi resolution 
classical grids. A multi-resolution classical grid at resolution level r has 2™ 
points. Moreover, it contains all the points of all the resolution levels r' < r. 
Every grid point Pj, in a multi resolution classical grid at resolution level r, can 
be written as Pi = . . . , |?} with < a\, . . . , a n < 2 r — 1, <zi, . . . ,a n G N, 
with the corresponding grid region (that we call as cells) £,(r,i) — [ ai '°^ +1 ) x 

[ a "-°; i+1 ), Here i is the unique identifier for each cell £, and it denotes its 
order in the generated samples. For any resolution level j, the resolution level 
j + 1 satisfies £j+i = £j/2. For any two cells £1,^2, at resolution levels ri,r 2 
respectively, with r\ < r 2 , and, £1 n £ 2 7^ 0, £1 will be called as the parent of £ 2 
at resolution ri, and, £ 2 as a c/iiZd of £1 at resolution r 2 . 

Ideally, one would like to use an ordering of samples that minimizes the 
discrepancy to find a counter-example as soon as possible. But discrepancy- 
optimal orderings take exponential time to compute, and exponential space to 
be stored (in n; see [23]). For our work, we use orderings that maximize the 
mutual distance between sampled grid points (see [23]). The mutual distance 
of a set K, is defined as p m (/C) = min^^g^; p(x, y). These orderings can be 



represented by using generator matrices, that are binary matrices of size n x n 
and represent bijective linear transformation over Z2. Any sample with identifier 
i, at resolution j, can be generated by using the ordering recursively on the bit 
representation of i. This method can also be used to bias the search towards the 
unsafe set T, by possibly changing the generator matrix. In Fig. 4, we show the 
samples for few resolution levels in 2 dimensions. As can be seen, cells whose 
identifiers are close to each other (e.g. Resolution = 2 or 3, i = 0, 1) are spaced 
quite far apart. This is in sharp contrast to the ordering that one would get by 
using a naive sampling scheme, like scanning for example. 



(a) Resolution— 1 



(b) Rcsolution= 2 



(c) Resolution— 3 



Figure 4: Ordering of samples based on mutual distance 



6 Algorithms for Resolution-Complete Safety Fal- 
sification 

In this section, we first explain the procedure for incremental construction of 
trajectories used by the algorithms. We then present the algorithms, and in 
Section 6.4, prove their resolution completeness. 

6.1 Incremental construction of trajectories 

As stated in Section 4.1, the algorithms use piece- wise constant control func- 
tions satisfying Proposition 2. As discussed in the proofs of Propositions 2, 3, 
for a given discrete state q £ Q, the dynamics can be simulated by using 
x{k) = AqXo + v, where xq £ Z(q, ■), and, the set of feasible inputs is given 
by v = T q u, u £ R km , \\u\\ < lj(q) = 1 - e 3 (q) \\T q + 1| . Note, that to be less 
conservative, we are making the input bound lj(q) dependent on q £ Q and the 
discretization Sj(q). Hence, the algorithms use the input bounds based on the 
space discretization and dynamics to incrementally build trajectories k steps 
at a time, by formulating the fc-step reachability problem as a linear program. 
To solve these linear programs more efficiently, we use ideas from [26] for multi 
parametric linear programming. 

We use the following additional notation in the remaining of this section. 
G = VJ q( zQG(q) denotes the union of all location specific grids in the algorithm. 
jo represents the smallest j, such that lj (q) = 1 — £j (g)||r 9 + || > 0,V<? £ Q. 



For all j > jo, and, Vq G Q, lj(q) = 1 — £j(<7)||r g + ||. e.j is a Q-dimensional 
vector with each element denoting value for a given location q. Since lj(q) can 
be different for each discrete location q E Q, at termination of the algorithms, 
completeness will be guaranteed with respect to lj — min gS g lj(q)- 



6.2 Breadth First Search with Branch and Bound 

The proposed algorithm is shown in Fig. 6. We will refer to this algorithm as the 
BFS-BB Safety Falsification algorithm. A few iterations of the algorithm 
for a first order discrete-time LTI system with single discrete mode, are shown 
in Fig. 5 for a maximum resolution level of j = 3. Cells marked in yellow are 
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Figure 5: An execution of BFS-BB Safety Falsification algorithm 



the ones that need to be explored further based on branch and bound strategy. 
Cells marked in red are conclusively not reachable from (where £ Gf is S 
at the left and £(r = 2,i = 0) at the right). Cells marked in blue are the ones 
that are found to be reachable. 

BFS-BB Safety Falsification (H , j max , G) 

1 {locations. Found} <— {0, 0} 

2 {jo , Gf , locations} <— init(-#.*S) //Initialization stop 

3 J v jo 

4 while (j < j max A — 'Found) do 



5 for all (<? G locations) do 

6 update _explorer (.q,j) //Update simulation parameters for current j,q 

7 for all (£ f e Gf(q) A ^Found) do 

8 for all (r G {1, 2, . . . J}) do 

9 for all (i e {0, . . . , 2 rT1 - 1}) do 

10 if (Found) then 

11 return (G/, Found, €j) / /Terminate if unsafe 

12 £t ■*— generate_sample(i, r, g) //Sample generation 

13 if (check_feasible(£t)) then 

14 success < — expand(£^ , £ t ) //Expansion step 

15 if (success) then 

16 add_node(£/, £t) //Update of G f (q) 

17 Found < — check_unsaf ety(^t ) //Check for unsafcty 

18 £f .attempted <— true //Mark as attempted 

19 locations . erase (q) 

20 j + 1 < — G.ref ine(j) / /Refinement step 



21 return (G / , Found , e j) 



Figure 6: BFS-BB Safety Falsification algorithm 



Data structure: Each cell £ has its identifier i, resolution level r, and 
boolean variables attempted , filled . attempted is true, if, £ 6 G/(q) has been 
attempted for expansion (and the cell is called attempted). The variable filled 
is true (and the cell is called as filled) if, the cell £ has size (j(q) (for some j) 
and is found to be reachable using lj(q), or, all its children at resolution j + 1 
arc filled. G/ is implemented as a /ias/i map to enable quick look up for existing 
cells in G/. 

Initialization step: In the first step, the algorithm computes the first 
feasible resolution level j and initializes Gf(q) Vq G Q. Next all q G Q are 
added to locations for which G/ 7^ 0. This happens in ±n±t(H.S) function in 
the algorithm. 

Exploration step: The grid is searched for reachable cells with current 
bounds on input in a recursive Breadth-First-Search (BFS) fashion along with 
Branch and Bound (BB) strategy in the algorithm for each location q G locations 
(after updating the simulation parameters used by the algorithm for location 
q in update_explorer(q, j) function in the algorithm). First, a filled cell £/ 
is chosen in G/(g). Next, for all depths r G {1, . . . , j} samples £t are gener- 
ated in generate_sample(i, r, q) function (after checking if Found = 0). The 
function check_f easible(£ t ) checks feasibility of £t based on branch and bound 
condition from coarse resolution levels, and the fact that its parent might have 
been marked as filled already at some coarser resolution level r' < r. The 
expand(£/, £ t ) function attempts to solve the /c-step reachability problem as 
discussed in Section 6.1. If it is successful then the add(£/,£ t ) function adds £/ 
to G f (q). 

Unsafety check: Each newly added cell £ t in add(£y , £ t ) step is checked for 
intersection with the unsafe set T, if, it is filled, in check_unsaf ety(£ t ) function. 

Discrete transition step: For a given location q, each cell £ G G(q) that 
is found to be reachable from G/(g) is checked for all the outgoing discrete 
transitions from q. If a guard G(q,q') is found to be enabled, then Q{q,q > ) H £ 
is added as a filled cell to Gf(q') using the function add(£, Q(q, q') n £) and q' is 
added to locations. This happens internally in add(£y,£ t ) function. 

Refinement step: When all the cells £f £ Gf(q) have been explored for 
expansion Vo G Q, and locations = 0, then the grid resolution is changed using 
the relation Ej+i(q) =Ej(q)/2, and the input bounds are changed from lj(q) to 
lj + i(q). For all the cells £ G Gf(q), and, Vq G Q, the variable attempted is reset 
to false. This happens in the G.refine(j) function in the algorithm. 

Termination criteria: To have a finite termination time, the refinement 
procedure is allowed only till j < j ma x, and, no counter example has been found. 

6.3 Resolution-complete Co-RRT 

In this section, we discuss a modified version of the Co-RRT algorithm proposed 
by us in [13] for continuous systems that is resolution-complete. This algorithm 
can be used for analyzing safety of a continuous system when it starts from 
equilibrium under the effect of exogenous inputs. We first state an important 



lemma regarding reachability of points that are convex combinations of two 
reachable points. 

Proposition 4. For an origin-stable, continuous system H , if S — {0}, then for 
any two reachable points ^i(0, U\, fci), ^2(0, U2, k%), u i> u 2 G U,k\,k2 G N, ki > 
k\, the convex combination ip\ — Xif^i + (1 — A)*02i ^ G [0, 1] is also reachable in 
ki time steps. 

Proof. Vi(0,ui,fci) = [BAB ...A k ^- l B]u x , V> 2 (0,wi,fc 2 ) = [£M_B . . . A k2 ~ 1 B]u 2 
u u u 2 G R kim , R k27n respectively. Let B k = [BAB . . . A k - 1 B], and, 6 C e R c , de- 
note the zero vector. This implies that ip\ = XBk 2 [ u i#m(fc 2 _fc 1 )]'+ (1 — A)5fc 2 u 2 . 
Since U is convex, w A = ^[ u 'i^m(k2-ki)\' + ^ — A),w 2 eW. ■ 

Corollary 1. 1Z(U) is a convex set for a system H (as in Proposition 4)- 

The proposed algorithm, called as the Co-RC Safety Falsification al- 
gorithm, is shown in Fig. 8. In Fig. 7, we show a few iterations of the algorithm 
for a second order system. 
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Figure 7: An execution of Co-RC Safety Falsification algorithm 

Data structure: The data structure is the same as one in previous algo- 
rithm except that G/ now contains a special node G/.Hull that contains the 
information about the convex hull and which is used for actual expansion step, 
and unsafety checks. 

Initialization step: In the first step, the algorithm computes the jo, and, 
G/.Hull is initialized to S. This is done in init (H.S) function in the algorithm. 

Exploration step: The sample generation and search strategy are similar 
to previous algorithm. The function check_f easible(£ t ) checks feasibility of £ t 
based on the fact that £t may be contained inside G/.Hull, i.e, £ t c G/.Hull, 
or, its parent may not be reachable with the current G/. In such a case, 
the function check_f easible(^) returns false. If it returns true, then the 
expand(G/.Hull, £ t ) function attempts to solve the fc-step reachability problem 
as discussed in Section 6.1. If it is successful, then the Co(G/.Hull, £ t ) func- 
tion updates G/.Hull in the function Co(G/.Hull, £/). The operator Co updates 
G/.Hull to the convex combination of G/.Hull and vertices of the cell £j. 

Unsafety check: Each time a new cell is added to G/, intersection of 
G/.Hull and T is checked in the function check_unsaf ety(G/.Hull). 



Co-RC Safety Falsification (H , j max , G) 



1 Found <— 

2 {jo , G/ .Hull} < — G.iriit(_f/.5) //Initialization stop 

3 j < — jo, change < — false 

4 while (j < j max A — >Found) do 

5 update_explorer(j) //Update the simulation parameters for current j 

6 repeat 

7 change < — false 

8 for all (r (E {1, . . . do 

9 for all (i G {0, . . . , 2 rn - 1}) do 

10 if (Found) then 

11 return (G f , Found , c j) //Terminate if unsafe 

12 £t < — generate_sample(i, r) //Sample generation 

13 if (check_f easible(^t)) then 

14 success < — expand(G/ .Hull, £t ) //Expansion step 

15 if (success) then 

16 change < — true 

17 G/.Hull <— Co(G/.Hull, £ t ) //Update the Hull 

18 Found < — check_unsaf ety(G/ .Hull) // Unsafcty check 

19 until (^change) 

20 j + 1 < — G.refine(j) //Refinement step 

21 return (Gf , Found, €j) 



Figure 8: Co-RC Safety Falsification algorithm 

Refinement step: When all the cells at a given resolution level j have 
been explored (possibly repeatedly) for reachability from G/.Hull and no new 
cell can be reached, the grid resolution is increased according to the relation 
£j+i = max(£j/2, £ m in), and the input bounds are changed from lj to Z^+i.This 
happens in the G.refine(j) function in the algorithm. 

Termination criteria: To have a finite termination time, the refinement 
procedure is allowed only till j < j max , and, no counter example has been found. 

We next discuss the resolution completeness of the proposed algorithms. 

6.4 Resolution completeness 

We first prove two important lemmas required to prove the main result in The- 
orem 1. The first lemma proves that for a given maximum resolution level 
j > jo 3 , the algorithms terminate in finite time, and the second lemma proves 
the required set inclusion for resolution completeness. 

Lemma 2. Consider a LTI discrete-time hybrid system H and a given j max > 
jo- Then the Safety Falsification algorithms terminate infinite time. 

Proof. Consider the CoRC Safety Falsification algorithm first. The algo- 
rithm starts with the value Sj initially. If at any stage a counter example 
is found, finite time termination is trivially guaranteed. Otherwise, note that 
by Assumptions 1, 2, G is guaranteed to have a finite volume. Since there 
is a given bound jmax on the resolution, it is guaranteed that the algorithm 
will be able to refine the resolution only a finite number of times. The sum 
of the number of cells over all the possible refinements of the grid G is given 

cyn [eg / e j max 1 1 

by -/Vcdis = 2 n -T ■ Hence in the worst case, the algorithm terminates 



3 For j < jo the input bounds are infeasible by choice of eq as discussed in Section 6.1 



within N^ clls iterations. Now consider the BFS-BB Safety Falsification al- 
gorithm. We have \Q\ locations. For each location q, let -/V C eiis(ij) be the total 
number of cells over all possible refinements. The number of guards can be 
no more than \ Q\ 2 . Hence, in the worst case, the algorithm terminates within 
|Q| 2 max, eQ iV c 2 cns((?) iterations. 

Lemma 3. Consider a LTI discrete-time hybrid system H and a given j max > 
jo- Then the Safety Falsification algorithms find a feasible counter example or 
else generate an approximation 7lj max such that 7?°(Wj max ) C 1Zj mso C 1Z(U). 

Proof. For a given j max , lj max is fixed. The set inclusion 72.j max C follows 
from the discussion in Section 3, and Proposition 2, 3. This guarantees feasibil- 
ity of the counter examples found by the algorithms. For a location q G Q, if a 
point is found reachable by the algorithms in a cell £, then the algorithms mark 
that cell as reachable based on the relaxation of input bounds from lj(q) to 1, 
where jo < j • < j max - As a result we are guaranteed that when the algorithms 
terminate, the set inclusion 7?.°(Z/ ?max ) C 72.j max also holds true. Taking convex 
combinations of reachable cells in the CoRC Safety Falsification algorithm 
does not violate this inclusion. ■ 

Theorem 1. Consider a LTI discrete-time hybrid system H as in Lemmas 2, 3. 
Then the Safety Falsification algorithms are resolution complete for the system 
H. 

Proof. Let j ma x > jo be given. Proposition 1 guarantees the existence of re- 
quired sequence of family of control functions. From Lemma 3 we know that if 
a counter example is found it is feasible. If no counter example is found, then it 
implies that there doesn't exist one (using the class of control functions Wj max ), 
from the set inclusion proved in Lemma 3. Moreover,, Lemma 2 guarantees that 
the algorithms will terminate in a finite time. I 

7 Simulation Results 

In this section, we present the simulations results obtained on different discrete 
time systems. The implementation has been carried out in C++ on a Pentium 
4, 2.4 GHz machine, with 512 MB of RAM. We will examine the performance 
of algorithms presented in Section 6.2 (called as BFS-BB), Section 6.3 (called 
as CoRC) and the one presented in [21], which is based on Depth-First-Search 
with Branch and Bound strategy (called as DFS-BB). 

7.1 Safety falsification of a discrete-time fifth-order sys- 
tem 

The example presented in this section is mainly intended to investigate perfor- 
mance of different algorithms for problems in moderate dimensions. We con- 
sider a fifth-order system, with dynamics specified as: x(i + 1) = Ax(i) + Bu(i), 
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Figure 9: Performance comparison of different algorithms for fifth-order system 



A = 0.6065 I, B — 3.935 I. The initial set is S = {0}. The simulation parame- 
ters are as follows: e = 27.2, j max = 4. This corresponds to ||w|j < 0.5679. The 
unsafe set is a hypercube of size 0.90 and is given by T = x + [— 0.45, 0.45] 5 . 
We did 50 test runs each with an unsafe set randomly centered at x, where 
x G [—7, 7] 5 was chosen using a uniform distribution. 

The performance results using different algorithms are shown in Fig. 9. The 
performance results indicate that when ever a counter example was found, all the 
algorithms terminated within 5 minutes. The BFS-BB algorithm falsifies safety 
sooner than other two algorithms. From the results, it is difficult to comment if 
CoRC performs better than DFS-BB or not. We have also done some profiling 
of the CoRC algorithm, and have found that almost 70% of the time is spent in 
removing the redundant vertices of the hull in the Co(G/.Hull, £/) function. In 
our current implementation, we reconstruct the convex hull from scratch every 
time the Co(G/.Hull, £y) function is called. We are working on using software 
libraries that support incremental construction of convex hulls. 

7.2 Safety falsification of a discrete-time second-order hy- 
brid system 

We now consider a second-order hybrid system with two discrete states. The 
example is interesting because of the fact that the guards and the invariants 
are all the same, and hence there is a potential problem of cycles. The dy- 
namics are specified as x(i + i) — A q x(i) + B q u(i),q G {1,2}, with, A\ — 
[0.679 0.404; -0.674 0.140], B x = [0.440; -0.213], A 2 = [0.679 -0.404; 0.674 0.140], 
B 2 = [0.3486; -0.1628]. T(l,-) = 1(2, ■) = 5(1,2) = 0(2,1) = [-2,2] 2 . 
The initial set is S = x [—0.1248, 0.1248] 2 . The simulation parameters are 
£o = 4, j max = 7. Corresponding to this ||tt|| < 0.792 for q — 0, and ||u|| < 0.869 
for q = 1. Hence, for the case when no counter example is found, completeness 
will be guaranteed for ||u|| < 0.792. 

The simulations are run for two cases based on the size of the unsafe set. 
In Case I, the unsafe set is small which makes it more difficult to be found 
by the explorer. However, quite often the unsafe set is specified as a large set 
(e.g. a half space representing the separation between two cars or aircrafts). To 




Figure 10: Safety falsification example and performance for Case I 
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Figure 11: Safety falsification example and performance for Case II 



evaluate the performance of the algorithm for such cases, we present test runs 
in Case II. We did 200 runs for both the cases. 

Case I: Small unsafe sets: For this case, the unsafe set is given by: T = 
q r x {x + [— 0.1, 0.1] 2 }, where x € [— 1.5,1. 5] 2 , was chosen using a uniform 
distribution and q r € {0, 1}. 

Case II: Large unsafe sets: For this case, we chose the unsafe set as a 
randomly oriented half space, given by: T = q r x {[xx, X2) ■ mx\ + ax-i < c}, 
where, m £ [—3.73, 3.73], c e [— 1.9,0], q r £ {0,1}, and a — — sgn(c). c, m are 
chosen from a uniform distribution such that 5flT = 0. 

Fig. 10 and Fig. 11 show performance results for the two cases along with 
a run when a counter example was found. Both the DFS-BB and BFS-BB 
algorithms falsify safety sooner for Case II. The BFS-BB falsifies safety in less 
than 10 seconds for 68 cases, compared to 54 for DFS-BB in Case I, and 164 
times compared to 149 times for Case II. 



8 Conclusions 

In this paper, we have presented sampling-based resolution-complete algorithms 
for safety falsification of linear time invariant discrete-time systems over infinite 
time horizon. The algorithms attempt to generate a legitimate counter example 
by incrementally building feasible trajectories in the state space at increasing 



levels of resolution or provides a guarantee in a finite time that no such example 
exists, when the input is restricted to a certain class. As an additional result, 
when no counter example is found, the algorithms provide us with an arbitrar- 
ily good under approximation to the reachable set whose quality is independent 
of length of trajectories. Efforts are currently underway to develop more ef- 
ficient algorithmsthat combine the nice features of both the depth-first-search 
and the breadth-first-search strategies to explore the state space. We are also 
investigating if its possible to develop similar algorithms for nonlinear systems. 
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